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Determining past changes in the amplitude, frequency and teleconnections of the El 
Nino-Southern Oscillation (ENSO) is important for understanding its potential sensitiv- 
ity to future anthropogenic climate change. Palaeo-reconstructions from proxy records 
s provide long-term information of ENSO interactions with the background climatic state 
through time. However, it remains unclear how ENSO characteristics have changed 
through time, and precisely which signals proxies record. Proxy interpretations are un- 
derpinned by the assumption of stationarity in relationships between local and remote 
climates, and often utilise archives from single locations located in the Pacific Ocean to 
io reconstruct ENSO histories. Here, we investigate the stationarity of ENSO teleconnec- 
tions using the Last Millennium experiment of CMIP5 (Coupled Model Intercomparison 
Project phase 5) (Taylor et al., 2012). We show that modelled ENSO characteristics 
vary on decadal- to centennial-scales, resulting from internal variability and external 
forcings, such as tropical volcanic eruptions. Furthermore, the relationship between 
is ENSO conditions and local climates across the Pacific basin varies throughout the 
Last Millennium. Results show the stability of teleconnections is regionally dependent 
and proxies may reveal complex changes in teleconnected patterns, rather than large- 
scale changes in base ENSO characteristics. As such, proxy insights into ENSO likely 
require evidence to be synthesised over large spatial areas in order to deconvolve 
20 changes occurring in the NIN03.4 region from those pertaining to proxy-relevant local 
climatic variables. To obtain robust histories of the ENSO and its remote impacts, we 
recommend interpretations of proxy records should be considered in conjunction with 
palaeo-reconstructions from within the Central Pacific. 
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1 Introduction 

The El Nino-Southern Oscillation (ENSO) is an important determinant of climate vari- 
ability, altering global rainfall patterns and modulating global temperatures. Under- 
standing the long-term characteristics of ENSO variability and its sensitivity to external 
5 forcings, such as greenhouse gases, represents a fundamental climate modelling and 
data challenge. Changes in ENSO behaviour may occur under future global warming 
(Power et al., 2013). However, there is a is a large dispersion in global climate model 
(GCM) projections of changes in ENSO characteristics (e.g. Collins et al., 2010; Vecchi 
and Wittenberg, 2010), and hence the sensitivity of the coupled ocean-atmosphere 
io system to future changing boundary conditions remains uncertain (DiNezio et al., 

2012) . In addition, investigations of the sensitivity of ENSO to anthropogenic climate 
change are restricted by the relatively short instrumental record, which provides us with 
limited guidance for understanding the range of ENSO behaviours. For example, the 
observed changes in the character of ENSO since the mid-1 970s towards a dominance 

is of El Nino, rather than La Nina, episodes are difficult to evaluate in terms of a forced 
response or unforced variability given the limited observational record almost certainly 
does capture the full range of internal climate dynamics. 

High resolution palaeo-reconstructions, including from tree rings, sediment cores, 
corals and speleothems, have the potential to provide long-term information about 
20 changes in modes of climatic variability and their sensitivity to different boundary con- 
ditions. Some tropical proxy records reveal ENSO interactions with the background 
mean climatic state. For example, data from long-lived fossil corals are often inter- 
preted quantitatively as estimates of ENSO changes through time that show a range 
of ENSO frequencies and amplitudes through time. Central Pacific coral reconstruc- 
25 tions generally reveal a weakened ENSO during the early Holocene (McGregor et al., 

2013) and highly variable ENSO activity throughout the Holocene (Cobb et al., 2013), 
which may have arisen from internal ocean-atmosphere variability (Cobb et al., 2003). 
Developing robust estimates of natural ENSO variability over a period longer than per- 

1581 


mitted through the instrumental record is a useful research avenue, with the potential 
for informing meaningful adaptive strategies for future climate change. 

Palaeo-ENSO proxy records of the Last Millennium (1000 years) are sparsely popu- 
lated temporally and spatially, and reconstructions remain uncertain (Cobb et al., 2003; 
5 Khider et al., 2011). It also remains unclear as to precisely which climatic signals as- 
sociated with ENSO are being recorded in these individual proxy records and whether 
these provide the necessary resolution to reconstruct ENSO changes. The assump- 
tion of stationarity of relationships between local and remote climates (teleconnections) 
underpins the interpretation of many palaeoclimate reconstructions, although station- 
io arity should not necessarily be assumed in terms of ENSO variability (Gallant et al., 

2013) . Are palaeo-reconstructions from the tropical Pacific recording base changes 
in the ENSO system or rather changes in teleconnected patterns? Previous model- 
based studies have identified sensitivity in the relationship between ENSO and the 
background climate state, and urged caution in the reconstruction of ENSO from proxy 

15 records under the assumption of stationarity of observed teleconnections (Coats et al., 
2013; Gallant et al., 2013). 

However, these studies have not comprehensively addressed the degree to which 
uncertainty about the non-stationarity of ENSO teleconnections can be assessed for 
particular locations and for particular mean climatic states. Furthermore, although we 
20 previously investigated the potential non-stationarity of hydrologic responses to ENSO- 
like conditions under disparate boundary conditions in idealised model simulations, we 
did not provide specific guidance for interpreting tropical proxy records (Lewis et al., 

201 4) , which currently comprise our dominant source of information about ENSO char- 
acteristics beyond the instrumental record. In addition, while previous studies have 

25 utilised proxy records, together with simulations using global climate models (GCMs) 
to evaluate the representation of ENSO in the current generation of GCMs (Cobb et al., 
2013), these approaches focused on using palaeo-ENSO reconstructions to test the 
performance of GCMs for the purpose of constraining uncertainty in future projections 
of ENSO behaviour under climate change. 
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As such, precisely which expressions of ENSO are being recorded in proxy archives 
under differing climatic boundary conditions have not been comprehensively interro- 
gated. Climate models, in addition to observational and proxy climate evidence, allow 
an understanding of long-term ENSO changes through time to be obtained (Schmidt, 
5 2010). A new generation of climate models and experiments has recently become 
available (Taylor et al., 2012), providing an opportunity for the first time to investi- 
gate ~ 1200 years of ENSO variability and establish a framework for understanding 
ENSO changes through time, using more models than previously possible. Hence in 
this current study, we investigate changes in ENSO characteristics (frequency and am- 
io plitude) in model experiments of the Last Millennium (“pastlOOO”). Focusing on three 
key climatic regions (East, Central and West Pacific), where explicit palaeo-ENSO re- 
constructions have been made, teleconnected patterns (the relationship between local 
and remote climates) throughout the Last Millennium are examined for surface tem- 
peratures and precipitation. We ultimately aim to determine whether proxy archives in 
is the tropical Pacific are likely to be recording alterations in ENSO base frequencies or 
local-scale teleconnections under differing boundary conditions. 


2 Datasets and methods 
2.1 Definitions 

The study is primarily focused on palaeo-ENSO variability from the tropical Pacific. 
20 Model data were investigated in three regions that have been identified as sensitive to 
modern ENSO variability and have also been used explicitly to reconstruct past ENSO 
changes (e.g. Cobb et al., 2013; McGregor et al., 2013). Area-mean anomalies for 
precipitation and surface temperature were calculated for the West (10° S-10° N, 105- 
1 55° E), Central (1 0° S-1 0° N, 1 70-1 30° W) and East Pacific (20° S-5° N, 65-90° W) re- 
25 gion and surface temperature for the NIN03.4 region (5° N-5° S, 170-120° W) (Fig. 1). 
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These regions are not intended to provide exhaustive coverage of ENSO impacts, but 
are large enough to provide useful comparisons with model-based data. 

El Nino episodes were defined based on surface temperature anomalies in the 
NIN03.4 region, with events defined in the models when NIN03.4 temperature anoma- 
5 lies were > 0.5 K for at least six consecutive months (Trenberth, 1 997). Conversely, La 
Nina episodes were defined when NIN03.4 temperature anomalies were < -0.5 K for 
at least six consecutive months. Spatial patterns are examined by compositing temper- 
ature and rainfall anomalies into positive (El Nino) and negative (La Nina) phases using 
these definitions for all CMIP5 models analysed. We utilise the NIN03.4 region as an 
10 index to classify ENSO conditions. Although the NIN03.4 region is commonly used to 
categorised ENSO episodes, it should be noted that there are other indices of ENSO 
that may also provide useful information beyond the Central tropical Pacific conditions 
described by the NIN03.4. 

2.2 Models and experiments 

15 CMIP5 data (Taylor et al., 2012) were downloaded from the Project for Model Diagno- 
sis and Intercomparison (PCMDI) through the Earth System Grid (ESG). Simulations 
were used of the historical (1 850-2005 CE) experiment, which is forced using chang- 
ing atmospheric compositions due to observed anthropogenic and volcanic influences, 
solar forcings and emissions of short-lived species from natural and anthropogenic 
20 aerosols. In addition, simulations were used of the Last Millennium (past 1000) (850- 
1850CE), in which reconstructed time-evolving exogenous forcings are imposed, in- 
cluding changes in volcanic aerosols, well-mixed greenhouse gases, land use, orbital 
parameters and solar changes. Each model’s pre-industrial control simulation (piCon- 
trol) with non-evolving pre-industrial forcings was analysed. Models were analysed 
25 where past 1000 simulations were provided (Table 1). One model (bcc-csml-1) was 
excluded from analysis because its representation of ENSO spectra is too short, pro- 
ducing a biennial oscillation rather than a broad spectral peak in the 2-8 year band. 
The MIROC-ESM model was also excluded, as it exhibits large drift related error in the 
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form of long-term trends that cannot be attributed to natural variability (Gupta et al., 
2013). Further details are provided in Sect. 3 on Model Evaluation. 

Data (precipitation (pr) and surface temperature (ts)) for six remaining models were 
regridded onto a common 1 .5° latitude by 1 .5° longitude grid. For the piControl and past 
s 1000 experiments, monthly anomalies were calculated by subtracting the mean sea- 
sonal cycle for each model. For the historical experiment, only model years 1 976-2005 
are considered, as this experiment necessarily produces a non-stationary climate due 
to the time-evolving anthropogenic greenhouse gas forcings imposed. Additional exper- 
iments were analysed for CMIP5-participating models, where available. For GISS-E2- 
io R (Schmidt et al., 2014) and IPSL-CM5A-LR (Dufresne et al., 2013) models, extended 
control simulations of > 500 years in duration were analysed and compared to forced, 
pastlOOO experiments. 

2.3 Diagnosing ENSO changes and teleconnections 

The location of ENSO activity in the historical and Last Millennium experiments in was 
is first explored using the leading empirical orthogonal function (EOF) of the tropical Pa- 
cific surface temperature. These spatial patterns were compared to the NIN03.4 index 
to determine possible non-stationarities in the site of ENSO activity through time (Li 
et al., 2011). This EOF analysis (Fig. SI) demonstrates that in both experiments, the 
surface temperature patterns are loaded in the NIN03.4 region. Although there are 
20 some differences in the spatial patterns of the leading EOF mode across the equato- 
rial Pacific, the similarity in model experiments in this particular region indicates that 
areal-average NIN03.4 temperatures provide a useful metric of ENSO activity in both 
experiments. 

Hence, a wavelet analysis was next used to examine the frequency and amplitude 
25 of NIN03.4 surface temperature variability in each model for statistically significant 
changes. Wavelet analysis is useful for examining non-stationary signal and provides 
time and frequency localisation. A Morlet mother wavelet (Torrence and Compo, 1998) 
with degree 6 was used to calculate the wavelet power spectra and identify large- 
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scale changes in variance. Wavelet spectral estimates were tested against red noise, 
represented as a first order autoregressive process. First, the wavelet power spectrum 
of NIN03.4 surface temperatures was considered, with variance in the wavelet power 
spectrum plotted as function of time and period (Fig. 2a). Next, the NIN03.4 mean 
s wavelet power spectrum, generated using a Morlet wavelet of 6°, was used as a metric 
for ENSO amplitude. The spectral power was calculated for the historical simulation 
and compared to the range of spectral power displayed in the pastlOOO experiment, 
calculated using ten 100 year epochs (Fig. 2b). 

The relationship between ENSO variability and teleconnected patterns in the tropical 
io Pacific regions (East, Central and West) was diagnosed through several complemen- 
tary approaches. First, an ordinary least squares regression between monthly NIN03.4 
mean surface temperature and remote area-mean surface temperature, and between 
monthly NIN03.4 mean surface temperature and remote area-mean precipitation was 
compared for the historical and Last Millennium experiments, for each region. Second, 
15 the relationship between local and NIN03.4 climates was considered using the corre- 
lation between variables (Corr(Local, Remote)), analogous to considering land-surface 
coupling strength (Lorenz et al., 2012). Correlations coefficients were calculated for 
monthly timeseries in ten 100 year epochs comprising the Last Millennium. Values 
were determined at each model gridbox and an area-weighted mean calculated for 
20 each region. The significance of correlations was assessed at the 95% confidence 
level for each coefficient using a t test. Third, the significance of identified changes in 
local-remote relationships during the Last Millennium was investigated. 

For each 100 year epoch comprising the Last Millennium, the El Nino- and La Nina- 
associated local temperature and precipitation anomalies were selected for each re- 
25 gion. A two-sided Kolmogorov-Smirnov (KS) test was used to investigate whether the 
distribution of local climate variables in 100 year epochs within the Last Millennium 
could statistically have been drawn from the same population (at the 5% significance 
level). A two-sided KS test was applied to each ENSO phase for each variable (sur- 
face temperature, precipitation) in each region (East, Central, West) comparing every 
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permutation of epochs sequentially (e.g. comparing El Nino-associated Central Pacific 
temperatures during 850-949 with 950-1049, then 1050-1149, then 1150-1249 etc.). 
A KS test was used for detecting changes in ENSO-remote climate relationships in Last 
Millennium timeseries as it is non-parametric and requires no assumptions to be made 
5 regarding the distribution of the data. A change is detected where the null hypothesis 
(that the distributions considered were drawn from the same population) is rejected at 
the 5% significance level. 


3 Model evaluation 

The basic properties of El Nino-Southern Oscillation (ENSO) simulated in Coupled 
io Model Intercomparison Project phase 5 (CMIP5) models (Taylor et al. , 2012) have been 
comprehensively examined, relative to observations, in previous studies (Bellenger 
et al., 2013; Guilyardi et al., 2012). Here, ENSO was examined through 6 metrics -1) 
ENSO amplitude (Nino3 sea surface temperature (SST) SD), (2) structure (Nino3 vs. 
Nino4 amplitude), (3) frequency (root mean square error of Nino3 SST anomaly spec- 
i5 tra), (4) heating source (Nino4 precipitation SD), (5) the amplitude of the ENSO bien- 
nial component (the ratio of the Nino3 SST anomaly timeseries power in the 3-8 years 
and 1-3 years bands) and (6) seasonality of ENSO (ratio between winter November- 
January over spring March-May average Nino3 SST anomalies SDs). Notably, some 
models have a tendency to represent ENSO with a too short period of about 2 years 
20 (bcc-csml-1, Fig. S2) or to have a spectral peak displaced towards longer periods of 
between 3 and 8 years (CSIRO-Mk-3.5 and MIROC5). Models that have contributed 
Last Millennium simulations to CMIP5, but are identified as erroneously biennial in 
terms of ENSO spectral power are excluded from this analysis. In addition, long-term 
drift in long model simulations must be considered. The MIROC-ESM model is ex- 
25 eluded from this analysis as it exhibits large drift related error in the form of long-term 
trends that cannot be attributed to natural variability, but instead relate to deficiencies 
in model physics and numerics (Gupta et al., 2013) (Fig. 3). 
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The remaining six models that contributed CMIP5 Last Millennium simulations (Ta- 
ble 1) were compared to twentieth century reanalysis data (20CR) (Compo and 
Whitaker, 2011) as a proxy for the observed climate that has been widely used to 
investigate observed ENSO-rainfall relationships (King et al., 2015; Klingaman and 
5 Woolnough, 2013). The CMIP5 models were compared to reanalysis precipitation and 
surface temperature for several ENSO-related characteristics over the common period, 
derived from the CMIP5 “historical” simulation. First, to investigate the model represen- 
tation of ENSO spatial patterns, the first empirical orthogonal function of the tropical 
Pacific surface temperature anomalies was calculated for 20CR reanalysis and CMIP5 
10 multi-model mean (MMM) EOF (Fig. 4a and b) and similarly for precipitation anoma- 
lies (Fig. 4c and d). This demonstrated the surface temperature and precipitation pat- 
terns are qualitatively similar for reanalysis and models; surface temperature patterns 
are generally of the same sign, although the meridional width of tropical temperature 
anomalies is narrower than in the reanalysis estimates, and simulated precipitation 
15 patterns are similar to the reanalysis estimate in the Central Pacific, although posi- 
tive anomalies are located too far westward in the CMIP5 MMM, compared with ob- 
served. In addition, the relationship between NIN03.4 surface temperature anomalies 
and global precipitation fields in reanalysis was compared to the CMIP5 MMM (Fig. 4e 
and f). The correlation coefficients between NIN03.4 temperature anomalies and local 
20 precipitation are generally of the same sign in simulated and reanalysis fields, includ- 
ing positive correlations in the Central and East Pacific and negative correlations in 
the west Pacific. These reanalysis-model comparisons are broadly insightful about the 
model representations of ENSO. 

It should be noted, however, that EOF analysis does not necessarily reveal modes 
25 that can be readily interpreted physically. In addition, only the MMM patterns are 
shown, and variations exist in the representations of ENSO in individual models. 
Nonetheless, key ENSO characteristics in the models, determined here and in previ- 
ous studies are physically plausible and hence useful for understanding potential non- 
stationarity of teleconnections. Furthermore, as an identical set of models is utilised 
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for each experiment, possible biases in ENSO representations in CMIP5 models are 
not considered prohibitive to investigating changes in the stability of teleconnections 
through time. 


4 ENSO during the Last Millennium 
s 4.1 ENSO characteristics 

Models demonstrate significant variance in the ENSO-relevant band (2-8 years) 
throughout the Last Millennium simulation (Fig. 2a). Several models also display 
significant low-frequency variability beyond the 8 year threshold; for example, the 
FGOALS-s2, GISS-E2-R and MRI-CGC3 models exhibit multi-decadal scale variability 
io in NIN03.4 temperatures during the Last Millennium simulations. These decadal- to 
centennial-scale El Nino- and La Nina-like episodes during the Last Millennium sim- 
ulations are evident in all models analysed here (Fig. 3). Models also demonstrate 
a range of ENSO amplitudes (Fig. 2b), with low spectral power shown, for example, in 
the GISS-E2-R and MRI-CGMC3 (Stevenson, 2012) models. There are differences 
is of ENSO amplitude between the historical and Last Millennium experiments within 
individual models. In the historical, the amplitude is generally weaker at all periods 
for MRI-CGMC3 and GISS-E2-R. The amplitude of higher ENSO-relevant periods (5- 
8 years) in the historical is largely outside the range exhibited in the Last Millennium 
(Fig. 2b). Previous model-based studies have also revealed strong inter-decadal to 
20 inter-centennial modulation of ENSO behaviour (Coats et al., 2013; Wittenberg, 2009) 
and warn that such modulation may not be fully revealed by the comparatively short 
instrumental climate record available and hence could attach large uncertainties to 
ENSO metrics diagnosed from short records. 

Low frequency modulation of Last Millennium ENSO characteristics may result from 
25 internal variability, or they may be externally driven. External forcings, specifically large 
tropical volcanic eruptions occurring between 1250 and 1600CE (Fig. S3), may pro- 
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duce decadal- to centennial-scale ENSO responses. Across the model ensemble in- 
vestigated here, La Nina-like post-volcanic responses appear to be simulated on multi- 
decadal timescales (Fig. 3). However, the response of ENSO to explosive volcanism is 
uncertain and may be dependent on the strength of the forcing (i.e. the size of the erup- 
s tion). On seasonal to annual timescales, model evidence suggests the radiative forcing 
due to volcanic stratospheric aerosols induces a La Nina episode that is followed by 
an El Nino episode after the peak of the forcing (McGregor and Timmermann, 2011). 
Proxy reconstructions derived from tree rings across the Pacific reveal similar ENSO re- 
sponses to those simulated, with anomalous cooling reconstructed in the East-Central 
io tropical Pacific in the year of volcanic eruption, followed by anomalous warming oc- 
curring one year after (Li et al., 2013). This pattern of response is likely sensitive to 
the specific set of eruptions investigated (Wahl et al., 2014). Previously, the associa- 
tion of eruptions and subsequent El Nino episodes has been demonstrated for forc- 
ings larger than that observed during the historical period for Mt Pinatubo (Emile-Geay 
15 et al., 2008). By comparison, when considering the strongest eruption of the Last Mil- 
lennium in 1258CE, the eruption likely produced a moderate or strong El Nino episode 
that interrupted the prevailing La Nina-like conditions induced by increased solar ac- 
tivity during the Medieval Climate Anomaly (Fig. S4). Low solar activity during the Little 
Ice Age (1450-1 850 CE) has been linked to generally cooler conditions, in contrast to 
20 a possible solar forcing of warmth during the Medieval Climate Anomaly (950-1 250 CE) 
(Schurer et al., 2014). Analysis of the frequency characteristics of various solar forc- 
ings imposed on the models during the Last Millennium simulations (Fig. S5) shows 
some decadal- and centennial-scale power, but this is not qualitatively consistent with 
variance in simulated ENSO conditions (Fig. 2). 

25 Alternatively, decadal- to centennial-scale modulation of ENSO behaviour may re- 
sult from internal ocean-atmosphere dynamics (Li et al., 2013), rather than a response 
to exogenous forcings. Previous studies suggest low-frequency ENSO characteristics 
may not be produced in long control simulations (Landrum et al., 2013), which ex- 
clude external climate forcing. When the properties of ENSO simulated in the control 
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simulations where compared to Last Millennium CMIP5 experiments, the smaller sub- 
set of extended control simulations available and analysed here exhibit qualitatively 
similar variability to that shown in the externally-forced Last Millennium experiment 
(Figs. 5-6). For the GISS-E2-R (Schmidt et al., 2014) and IPSL-CM5A-LR (Dufresne 
5 et al., 2013) models, extended unforced pre-industrial control (“piControl”) simulations 
of > 500 years in duration were provided as part of CMIP5. 

The smoothed mean annual timeseries of NIN03.4 region temperature anomalies 
during the control simulations (Fig. 5) exhibit qualitatively similar variability to the forced 
Last Millennium simulations (Fig. 3), including multi-decadal to centennial- scale El 
io Nino- and La Nina-like phases. Using wavelet analysis (Torrence and Compo, 1998), 
the amplitude and frequency of ENSO variability is investigated throughout the piCon- 
trol simulations. The models reveal some differences in low-frequency variability be- 
yond the 2-8 ENSO-relevant band, when comparing the control (Fig. 6) and the Last 
Millennium (Fig. 2). In GISS-E2-R, for example, there is weaker spectral power ex- 
15 hibited on multi-decadal timescales in the control experiment, compared with the Last 
Millennium, although some low frequency modulation is apparent in both models due 
to internal variability. This limited suite of GMIP5 models hints that ENSO variability 
on timescales that are poorly captured in the short instrumental record results from 
internal climate variability, and hence occurs in both piControl and Last Millennium 
20 experiments. 

Overall, there may be important difference in ENSO characteristics in the Last Millen- 
nium experiment, compared to the historical and piControl simulations, relating to both 
external forcings from large tropical volcanic eruptions or internal variability revealed on 
longer than inter-annual timescales. Differences between the piControl and Last Millen- 
25 nium experiments may represent an artefact of comparing a limited number of model 
realisations, or rather difference may reveal forced responses resulting from some com- 
bination large tropical volcanic eruptions of the Last Millennium or a time varying solar 
cycle (Emile-Geay et al., 2008). As only two modelling groups contributed extended 
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control simulations, we have not attempted a formal detection of the significance of 
differences in ENSO characteristics between forced and unforced simulations. 

Rather, we suggest that these differences highlight the considerable natural variabil- 
ity in the ENSO system that cannot be determined entirely from centennial and shorter 
5 records (Wittenberg, 2009). Proxy reconstructions providing “snapshots” of climate his- 
tory throughout the Last Millennium are often too short to resolve unforced and forced 
changes and these brief snapshots may reflect sampling of specifics ENSO character- 
istics within a short window. Longer reconstructions of spliced coral records, however, 
suggest ENSO responds to both low-frequency internal variability and external forcings 
10 (Cobb et al., 2003; McGregor et al., 2013). 

4.2 ENSO impacts and teleconnections 

Models show broadly similar global impacts associated with NIN03.4 regional temper- 
ature anomalies in the Last Millennium and historical experiments (Fig. 1). Composited 
patterns of global surface air temperature anomalies associated with positive (El Nino) 
15 and negative (La Nina) ENSO phases derived from all analysed models are spatially 
coherent across the experiments. Also, the Last Millennium patterns largely resem- 
bling those in the historical experiment. Although ENSO surface temperature anoma- 
lies across the Pacific are qualitatively similar, anomalies associated with the historical 
period (1976-2005) are generally of greater magnitude, particularly at remote loca- 
20 tions outside the equatorial Pacific, including over North America and the south Pa- 
cific. These differences in magnitude between the Last Millennium and the historical 
may relate to the differing boundary conditions during the historical period associated 
with anthropogenic forcings, such as long-lived greenhouse gases, or simply from the 
greater diversity of ENSO episodes represented in the longer Last Millennium simula- 
25 tion. 

Hence, we next explore the temporal stability of teleconnected ENSO patterns and 
ENSO impacts in different domains during the Last Millennium. The relationship be- 
tween NIN03.4 regional temperature anomalies and the mean local climate is exam- 
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ined in each analysed Pacific region (East, Central, West). Overall, the linear relation- 
ship between local and NIN03.4 surface temperature for the East Pacific and Central 
Pacific regions is consistent in the Last Millennium and historical experiments (Fig. 7). 
However, when considering data from the analysed models for the West Pacific region, 
5 this linear relationship is different for these experiments. Also, there is a statistically sig- 
nificant difference between experiments in the relationship between local precipitation 
and NIN03.4 conditions in the East Pacific. 

Complexity in the relationships between local (surface temperature and precipita- 
tion) and remote (NIN03.4 surface temperature) conditions during the Last Millennium 
10 simulation is further considered using the correlation between variables (Corr(Local, 
Remote)). This approach is analogous to considering land-surface coupling strength 
(Lorenz et al., 2012). We diagnose temporal stability using this correlation in ten 
1 00 year epochs that comprise the Last Millennium (Fig. 8). The strength of the remote- 
local relationship varies temporally and is regionally and climate variable dependent, 
is In the West Pacific, particularly, this coupling is generally weak and varies significantly 
between epochs. For local precipitation, the strongest coupling occurs for the Central 
Pacific. 

We also investigate the significance of these identified Last Millennium changes in 
local-remote relationship. A Kolmogorov-Smirnov (KS) test was used to determine 
20 whether the distributions of El Nino- and La Nina- associated local temperature and 
precipitation anomalies in each region in 100 year Last Millennium epochs could sta- 
tistically have been drawn from the same population. There are detectable differences 
(at the 5% significance level) in the distribution of ENSO-associated local climate vari- 
ables in these 100 year epochs. West Pacific El Nino- and La Nina- associated temper- 
25 atures, for example, significantly vary in character through the Last Millennium. Tempo- 
ral changes in local ENSO fingerprints (Corr(Local, Remote)) of the Last Millennium, 
also likely result from external forcings and/or internal ocean-atmosphere dynamics. 
However, these same relationships were not explored in the extended control simu- 
lations because of the small number of contributions available from different models. 
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Regardless of attributing the source of ENSO variability in the Last Millennium simu- 
lations to external forcings or internal dynamics, it is evident in the model experiments 
that differing teleconnections may result at different points in time and may differ from 
present-day relationships. In addition, Last Millennium variability in ENSO-local climate 
5 relationships across sites the tropical Pacific suggests that global ENSO changes do 
not necessarily scale linearly to local scales and cannot be assumed to do so. 

4.3 ENSO under differing boundary conditions 

The CMIP5 archive also provides simulations of the mid-Holocene (midHolocene, circa 
6000 years ago) from multiple participating climate models. The mid-Holocene pro- 
10 vides a well-constrained target for model-based studies (Schmidt et al., 2004) with 
substantially larger time-evolving forcings than those imposed during the Last Mil- 
lennium, and this period has also been the target of palaeo-reconstructions. Hence, 
these simulations are also briefly investigated here, in addition to the information pro- 
vided by the Last Millennium experiment. Mid-Holocene simulations are run for at least 
15 100 years after reaching equilibrium and have changed orbital parameters and atmo- 
spheric concentrations of greenhouse gases imposed. Other boundary conditions such 
as aerosols, solar constant, vegetation and topography are prescribed as the same as 
in the pre-industrial control simulation. 

By way of context, Cobb et al. (2013) report that Central Pacific corals record highly 
20 variable ENSO activity through the Holocene, although no systematic trend in ENSO 
variance was demonstrated in this study. A complementary Central Pacific reconstruc- 
tion from Kiritimati Island suggests that ENSO variance was persistently reduced by 
79%, compared with today at this location about 4300 years ago (McGregor et al., 
2013). Central Pacific coral-based evidence of ENSO variability is substantially differ- 
25 ent from lower-resolution records from the eastern equatorial Pacific (Conroy et al., 
2008; e.g. Moy et al., 2002). Collectively, East Pacific records suggest a systematic de- 
crease in mid-Holocene ENSO variance. On the West Pacific side of the basin, corals 
from northern Papua New Guinea reveal a reduction in ENSO frequency and ampli- 
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tude over the period of 7. 6-5. 4 ka (thousand years ago) compared with today, and also 
identifies large and protracted El Nino events for 2.5-1 .7 ka (McGregor and Gagan, 
2004). These Mid-Holocene ENSO reconstructions do not necessarily provide con- 
tradictory information, but may instead reflect geographic complexities (Carre et al., 
5 2014; Cobb et al., 2013). However, as proxy-based reconstructions from each of these 
regions have been used to infer changes in the same coupled ocean-atmosphere sys- 
tem, we also examine teleconnected ENSO patterns under these significantly different 
boundary conditions that characterise the mid-Holocene. 

In this study, we consider the subset of participating CMIP5 models with contributions 
io of mid-Holocene simulations (MRI-CGCM3, IPSL-CM5A-LR, FGOALS-S2, CCSM4) 
and find a general reduction in spectral power across ENSO-relevant frequencies 
(Fig. S6a) that has also been reported in model experiments of this period conducted 
prior to the release of CMIP5 (Chiang et al., 2009). This reduced spectral power in 
the ENSO band can be considered a metric for reduced ENSO amplitude (Stevenson, 
is 2012). Previous model and proxy-based studies have also hinted at subdued ENSO 
activity in the mid-Holocene. For example, early studies using simple numerical mod- 
els of the coupled ocean-atmosphere system by Clement et al. (2000) demonstrate 
increasing ENSO variability throughout the Holocene in response to time varying or- 
bital forcings. The impact of mid-Holocene orbital changes on ENSO variability has not 
20 been demonstrated comprehensively from proxy records. However, various fossil coral 
reconstructions indicate that there may have been reductions in ENSO variability in the 
mid-Holocene (Cobb et al., 2013). 

In addition, when CMIP5 midHolocene model data are composited into positive (El 
Nino) and negative (La Nina) phases, the magnitude of simulated mid-Holocene spa- 
25 tial patterns of ENSO impacts (Fig. S6b and c) are subdued, relative to the historical. 
The relationship between NIN03.4 mean surface temperature anomalies and regional 
(East, Central, West Pacific) temperature and precipitation was examined (see Fig. S7 
for West Pacific). This indicates that the relationship between West Pacific surface tem- 
perature anomalies and corresponding NIN03.4 temperature anomalies differs from 
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the midHolocene and historical simulations. The frequency of high and low local sur- 
face temperature anomalies in the West Pacific during El Nino defined conditions is 
also reduced in the midHolocene experiment compared with the historical (Fig. S7). 
The West Pacific surface temperature anomaly distributions for both El Nino and La 
5 Nina phases are notably narrower for the midHolocene, relative to the historical, indi- 
cating a reduced magnitude of ENSO impacts in this region. Conversely, the precipita- 
tion distributions for the West Pacific are for the most part similar for each experiment. 
The NIN03.4 impacts on East and Central Pacific regional temperatures are broadly 
similar for the historical and mid-Holocene. 


10 5 Towards reconstructing robust ENSO histories 

Palaeoclimate studies typically utilise archives from single locations located in the Pa- 
cific Ocean to reconstruct generalised basin-scale histories of ENSO. However, we 
demonstrate that proxies in one location alone should not be considered regionally 
representative, or singularly insightful about robust ENSO reconstructions without ex- 
15 plicit examination of the stability of ENSO teleconnections. Using standardised simula- 
tions performed with multiple participating climate models, we have demonstrated that 
there is model evidence for decadal- to centennial-scale ENSO variability during the 
Last Millennium that may not be fully revealed by the comparatively short instrumental 
climate record. The existence of varying ENSO characteristics throughout the Last Mil- 
20 lennium is also supported by proxy-based climate reconstructions (Cobb et al., 2003), 
which show variable ENSO characteristics include changing frequency and amplitude 
compared to modern during the Last Millennium. 

Furthermore, we find ENSO-local climate relationships vary significantly throughout 
the Last Millennium in the various tropical Pacific regions. In addition to (1) the sensi- 
25 tivity of ENSO to the background climate state and (2) the decadal- to centennial-scale 
modulation of ENSO behaviour from internal variability and external forcings such as 
volcanic eruptions, this suggests that (3) changes in ENSO do not necessarily scale 
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linearly to local scale impacts. Models suggest it may be inherently difficult to decon- 
volve variability in the NIN03.4 region and local-scale, teleconnected climatic change 
in the West and East Pacific regions. The West Pacific Warm Pool is likely sensitive 
to subtle shifts in the western extent of the warm (cool) tongue characterising positive 
5 (El Nino) and negative (La Nina) episodes. Hence, subtle changes in the Pacific basin 
may impact this region through several ocean-atmosphere mechanisms. 

We argue that proxy insights into change and variability in ENSO system are likely to 
be most robust when evidence is be synthesised over large spatial areas. That is, only 
incomplete information about temporal changes in a large-scale climate system can 
10 be provided by considering changes at a singular location (i.e. a time series of a cli- 
matic variable). Rather, considering multi-dimensional information in the form of spatial 
patterns of change through time is likely to yield more robust insights in large-scale 
systems. This provides a framework for enhanced interpretations of the invaluable in- 
formation of palaeoclimatic change provided by proxy records. For example, combined 
15 evidence from the West and Central Pacific is more likely to reveal the potentially sub- 
tle changes in ENSO-associated spatial patterns of temperature and precipitation per- 
turbations across the Pacific. For remote regions outside the equatorial Pacific, the 
non-stationarity of ENSO teleconnections is likely to be more problematic. These sites 
should be considered in conjunction with palaeo-reconstructions from within the Cen- 
20 tral Pacific basin, the so-called “centre of action” of ENSO (Cobb et al., 2013). We do 
not expect that this recommendation is dependent on the NIN03.4 metric used to de- 
fine ENSO utilised here. Indeed, such recommended spatially integrated approaches 
have already been undertaken and provide valuable information over the recent past 
(e.g. Li et al., 2013). It is also worth noting that under boundary conditions significantly 
25 different from present, such as the mid-Holocene, the stability of ENSO teleconnections 
is likely to be more variable, and hence potential non-stationarities in local-remote re- 
lationships require explicit consideration in proxy interpretations. 

This study has used palaeoclimate simulations conducted using a suite of CMIP5- 
participating models with various forcing to investigate changes in ENSO and its 
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teleconnections under differing boundary conditions (the Last Millennium and mid- 
Holocene). Such palaeo-modelling approaches do not attempt to replace proxy re- 
constructions, but rather combining multiple approaches can provide enhanced inter- 
pretations of reconstruction of past climate guiding our understanding of the most con- 
5 sistent physical explanations for change (Schmidt, 2010). In this study, we considered 
a subset of CMIP5 models that contributed palaeo-simulations and these necessar- 
ily contain systematic biases in ENSO representations (Power et al., 2013). In their 
study focused on understanding ENSO responses to volcanic forcings, Emile-Geay 
et al. (2008) suggested further forcing/response insights could be provided by GCMs 
io with realistic ENSO cycles and asked whether the current generation of models were 
up to the task. Deficiencies in our theoretical knowledge of ENSO and the difficulties 
in representing physically realistic ENSO cycles in GCMs (Guilyardi et al., 2012) are 
a limit on providing robust quantitative understanding of forced and unforced changes 
in the ENSO system. Nonetheless, existing model simulations remain useful for ex- 
15 amining palaeoclimates, despite their biases. These reveal spatially and temporally 
complex changes in ENSO and its teleconnected patterns during the Last Millennium. 
These complex spatial and temporal changes should be considered when developing 
robust proxy interpretations and ENSO histories in order that these are most useful for 
constraining future ENSO behaviour under greenhouse forcings. 

20 This study has identified several avenues for further research. First, additional tar- 
geted experiments within a single climate model that well represents ENSO spatial 
dynamics, particularly on the western extent of the warm/cold tongue, would provide 
further insight into the apparent complexity of ENSO impacts through time. In particu- 
lar, explicit investigation of the sensitivity of changes in the location of maximum base 
25 ENSO activity across the Pacific to varying boundary conditions would provide use- 
ful information about the range of variability of this coupled system. Furthermore, our 
present study did not comprehensively investigate the relative influences on various 
external forcings (solar and volcanics) and internal variability on ENSO characteristic, 
which would provide useful information for comparison with proxy records. 
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More direct comparisons between model output and proxy reconstructions can be 
provided by using a larger set of hydrological diagnostics and by employing pseudo- 
proxy techniques. In the first approach, the relationships between local climate pa- 
rameters and circulation, beyond the instrumental period, are examined using global 
5 climate models. Incorporating water isotope tracers into global climate models as a di- 
agnostic tool brings output one step closer to direct model-proxy comparisons. On- 
going research using water isotope tracers in GISS ModelE-R in experiments of the 
historical and Last Millennium period will provide further location-specific information 
about water isotope proxies and ENSO variability. Future studies focusing on water 
10 isotope variability may provide useful insights into changing ENSO impacts across the 
Pacific as these hydrological diagnostics have been shown to integrate regionally co- 
herent climatic information, compared with precipitation, which is highly variable (Lewis 
et al., 2010). Furthermore, in the second approach, a simulated time series intended 
to mimic actual proxy records (“pseudo-proxy”) is generated from a climate model sim- 
i 5 ulation (Anchukaitis and Tierney, 2013). The pseudo-proxy approach can be used to 
interrogate the necessary proxy density required for producing skilful regional climate 
field reconstructions and provide guidance on interpretations of reconstructions from 
particular locations (Smerdon, 201 1 ; Wahl et al., 201 4). 


The Supplement related to this article is available online at 
20 doi : 1 0.51 94/cpd-1 1 -1 579-201 5-supplement. 
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Table 1. Details of CMIP5 experiments and models analysed. Further details can be found 
through the Program for Climate Model Diagnosis and Intercomparison (PCMDI). 


Experiment 

Major forcings 

Years Analysed 

Models 

historical 

Time-evolving anthropogenic 

(greenhouse gases, aerosols, 
ozone) and natural (solar, vol- 
canics) 

1 976-2005 CE 

CCSM4, FGOALS-S2, GISS-E2-R, HadCM2, 
IPSL-CM5A-LR, MRI-CGCM3 

pastlOOO 

Time-evolving greenhouse gases, 
solar, volcanics, land use and or- 
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Figure 1. Composited anomaly maps for surface temperature (K) for CMIP5 models (left, El 
Nino episodes; right, La Nin a episodes) for historical (a and b) and past 1000 (c and d) exper- 
iments. Rectangular boxes indicate the West, Central and East Pacific regions. 
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Figure 2. Wavelet analysis of simulation of the Last Millennium, showing (a) Morlet wavelet 
coefficients shown by colours ranging from white (weak) to deep red (strong), with the contours 
enclosing regions of greater than 95% significance level. Colourbars are omitted since units 
are arbitrary, (b) Global mean NIN03.4 power spectrum (K 2 unit -1 frequency, black), relative to 
the red-noise (AR(1)) benchmark (dashed). The historical simulation is shown in black and the 
5th-95th percentile range across the Last Millennium shown by purple envelope, calculated 
using 100 year epochs. Spectral power was calculated using a Morlet wavelet of degree 6. 
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Figure 3. Running annual-mean surface temperature anomalies (K) over the NIN03.4 region 
(5°N-5°S, 170-120° W) for Last Millennium simulation for each participating CMIP5 model, 
showing long-term drift in MIROC-ESM model. Red/blue shading highlights departures from 
each model’s long-term mean. Running means were calculated using a 240 month triangle 
smoother. The vertical grey lines at years 1 258, 1 453 and 1 600 mark the timing of large volcanic 
eruptions of the Last Millennium. 
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Figure 4. (a-d) Comparison of leading patterns (standardised, first EOFs) of monthly vari- 
ability in surface temperature and precipitation for reanalysis (left: (a), surface temperature; 
(b), precipitation), CMIP5 models (b), surface temperature; (d), precipitation). CMIP5 historical 
patterns are the multi-model mean (MMM) of the first EOF of each individual model, (e and 
f). Spatial correlation coefficients between NIN03.4 index and 20CR precipitation (e) and the 
CMIP5 MMM (f). Stippling indicates Spearman’s rank correlations significant at the 95% level. 
Rectangular boxes indicate the East, Central and West Pacific regions. 
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Arbitrary Year 

Figure 5. Running annual-mean surface temperature anomalies (K) over the NIN03.4 region 
(5°N-5°S, 170-120° W) for extended piControl simulations conducted with GISS-E2-R and 
IPSL-CM5A-LR models. Red/blue shading highlights departures from each model’s long-term 
mean. Running means were calculated using a 240 month triangle smoother. Control simula- 
tions are spun up to quasi-equilibrium and run for ideally >500 years, providing an arbitrary 
timeseries of model internal variability. 
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Figure 6. Wavelet analysis of simulation of extended piControl simulations conducted with 
GISS-E2-R and IPSL-CM5A-LR models, showing (a), Morlet wavelet coefficients shown by 
colours ranging from white (weak) to deep red (strong), with the contours enclosing regions of 
greater than 95% significance level. Colourbars are omitted since units are arbitrary, (b) Global 
mean NIN03.4 power spectrum (K 2 unit _1 frequency, black), relative to the red-noise (AR(1)) 
benchmark (dashed). 
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Figure 7. Scatter plots of mean surface temperature anomalies (K) in the NIN03.4 region 
against area-mean anomalies in surface temperature (left) and precipitation (right) for the East 
(a, d), Central (b, e) and West (c, f) Pacific. Data points are shown for the historical (grey) 
and pastlOOO (blue). Lines of best fit are calculated for each experiment using ordinary least 
squares regression. 
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Figure 8. Area-mean correlation coefficients (R) of NIN03.4 and (a) surface temperature, (b) 
precipitation for the East (black square), Central (red cross) and West (blue cross). Data points 
show correlation coefficients calculated for ten 1 00 year epochs comprising the Last Millennium 
simulation. 
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